Nonlinear effects in charge stabilized colloidal suspensions 
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Molecular Dynamics simulations are used to study the effective interactions in charged stabilized 
colloidal suspensions. For not too high macroion charges and sufficiently large screening, the concept 
of the potential of mean force is known to work well. In the present work, we focus on highly 
charged macroions in the limit of low salt concentrations. Within this regime, nonlinear corrections 
to the celebrated DLVO theory [B. Derjaguin and L. Landau, Acta Physicochem. USSR 14, 633 
(1941); E.J.W. Verwey and J.T.G. Overbeck, Theory of the Stability of Lyotropic Golloids (Elsevier, 
Amsterdam, 1948)] have to be considered. For non-bulklike systems, such as isolated pairs or triples 
of macroions, we show, that nonlinear effects can become relevant, which cannot be described 
by the charge renormalization concept [S. Alexander et al., J. Chem. Phys. 80, 5776 (1984)]. 
For an isolated pair of macroions, we find an almost perfect qualitative agreement between our 
simulation data and the primitive model. However, on a quantitative level, neither Debye-Hiickel 
theory nor the charge renormalization concept can be confirmed in detail. This seems mainly to 
be related to the fact, that for small ion concentrations, microionic layers can strongly overlap, 
whereas, simultaneously, excluded volume effects are less important. In the case of isolated triples, 
where we compare between coaxial and triangular geometries, we find attractive corrections to 
pairwise additivity in the limit of small macroion separations and salt concentrations. These triplet 
interactions arise if all three microionic layers around the macroions exhibit a significant overlap. In 
contrast to the case of two isolated colloids, the charge distribution around a macroion in a triple is 
found to be anisotropic. 



I. INTRODUCTION 

In order to simplify the description of colloidal sys- 
tems, one often tries to determine effective interactions 
between the colloidal particles, thus integrating out the 
solvent's degrees of freedom This is not a trivial 
task, because in general the effective interactions depend 
on the thermodynamic state of the system, and one is 
often confronted with the problem of thermodynamic in- 
consistencies . A problem that is of particular interest 
is that of effective interactions between charged colloids 
(macroions) in an electrolyte solution. On a mean-field 
level, such a system can be described by the Poisson- 
Boltzmann (PB) equation Q. The linearized version of 
this equation is the basis of the DLVO theory for charged 
colloids 3, and we will refer to it in the following as 
the Debye-Hiickel (DH) limit of the PB equation [^. 
In the DH limit, the problem can be solved analytically 
and yields a screened Coulomb potential for the effective 
interactions between the macroions. The characteristic 
range of this potential is given by the Debye length 
which is controled by solvent properties such as the salt 
concentration. 

The physical picture of the DH description is rather ap- 
pealing: Due to the charge of a macroion, a layer of thick- 
ness is formed around it, consisting of oppositely 
charged microions (counterions), leading to a screening of 
the bare Coulomb interaction. Although the DH descrip- 
tion is only valid for weakly charged macroions and if cor- 
relation effects between the microions in the electrolyte 
solution can be neglected, it is tempting to character- 
ize also the effective interactions between highly charged 



macroions by a potential of screened Coulomb form. In- 
deed, this is the idea of the famous concept of charge 
renormalization that has been put forward by Alexander 
et al. ^ . It is based on the observation that in the frame- 
work of the so-called cell model (see below) the numerical 
solution of the nonlinear PB equation can be fitted far- 
ther away from the macroions' boundaries by a screened 
Coulomb potential with a renormalized charge Z^s < Z 
(with Z the bare charge of a macroion). Trizac et al. 
recently extended the numerical recipe of Alexander et 
al. by providing an analytical scheme to calculate Zcs 
as well as the effective screening length and the effective 
salt concentration. 

As already mentioned, DH theory is based on the lin- 
earized PB equation, which implies pairwise additivity 
of interaction energies. On the other hand, for highly 
charged macroions, nonlinearities imply the occurrence 
of many-body interactions and thus pairwise additivity 
does not hold. The simplest system, in which many- 
body effects could be expected, consists of three iso- 
lated macroions in an electrolyte solution. Indeed, such 
a system has been studied in a recent experiment us- 
ing scanned optical line tweezers 0; 0]- ^^^^ work, 
charge-stabilized silica particles with a diameter of about 
1 /im suspended in water were considered. It was possible 
to measure three-body interactions directly, and it was 
found that, in agreement with numerical solutions of the 
nonlinear PB equation |10(] , three-body contributions 
to the total interaction energy are attractive. 

Also Molecular Dynamics (MD) computer simulations 
have been used to investigate systems of "isolated" 
macroion pairs and triples 0, 0, ITsL [T3 . ITsI l . In these 



2 



studies, charged colloids were investigated in the frame- 
work of the so-called primitive model. In this model, a 
system of macroions, counterions and salt ions is consid- 
ered without explicitly taking into account the uncharged 
part of the solvent. Based on the primitive model, Al- 
lahyarov and Lowcn found that DH theory works well 
for a system of two macroions and, in agreement 
with experiment ^ ^ and PB theory , that three- 
body contributions are attractive in the case of three 
macroions 0. These authors also studied a system 
of two macroions, in which uncharged solvent particles 
were added to the electrolyte solution j^^. An inter- 
esting finding of this work was that the neutral solvent 
leads to a renormalized charge, which is smaller than 
the bare charge of the macroions, similar to the concept 
proposed by Alexander et al. |a . In a different simula- 
tion study by Tehver et al. |l4| . the counterions were 
introduced via density distributions in the framework of 
a density functional theory. Surprisingly, in the case of 
three macroions, no evidence for many-body forces was 
found, and the forces could be well described by DH the- 
ory. 

In this work, MD simulations are presented that tie 
in with the previous simulation studies. We consider 
systems of two and three highly charged macroions in 
a primitive model solvent. In the two-particle case, we 
check to what extent DH theory describes the effective 
interactions; thereby, the effect of nonlinearities is quan- 
tified. This is done for different amounts of added salt, 
focussing on the limit of small salt concentrations. In a 
second step, we address the influence of nonlinear effects 
on three-body interactions. In order to study these ef- 
fects, a triple of macroions is considered in two different 
geometries by placing the macroions on an equililateral 
triangle or along a straight line. We check on whether 
the concept of charge renormalization can also be applied 
to isolated pairs or triples of particles. Furthermore, we 
ask for the validity of the mean-field description and how 
effective interactions develop from the two-particle case 
to the bulk. Our major concern is the influence of non- 
linearity, which can be seen for high macroion charges 
and low salt concentrations. We especially consider cases 
of overlapping or interacting Debye-layers in the case of 
non-bulklikc macroion configurations. 

Our paper is organized as follows: After briefly dis- 
cussing some results of DH theory and the concept of 
charge renormalization, we give an overview of the simu- 
lation details. In Sec. lIVI ^ we present our results for sys- 
tems that consist of a pair of macroions, and in Sec. lIVB 
systems with macroion triples are considered. Finally, we 
discuss the results and draw some conclusions. 



II. DH POTENTIAL, PB EQUATION, AND THE 
CONCEPT OF CHARGE RENORMALIZATION 



e is the elementary charge and Z the valency). They 
are immersed into a polar, structureless medium with di- 
electric constant e. This medium is characterized by the 
Bjerrum length Ab = /{AirekBT), i.e. the distance at 
which the electrostatic energy between two point charges 
equals the thermal energy k^T. 

In the DH limit, the interaction potential between 
two macroions, separated by a distance r, is given by 
a screened Coulomb (Yukawa) potential 0, 



u{r) = /cbTAb 



Z exp(Kf7/2) 



1 + Kcr/2 



exp(— Kr) 



where 



K = \/47rAB(2ns + Znc)/V 



(1) 



(2) 



is the screening parameter, ric represents the number of 
macroions, and Ug is the number of added salt ion pairs. 
In Eq. 101, it is assumed that the electrolyte solution is 
formed by monovalent microions in a system of total vol- 
ume V. The microions consist of Zn^ negatively charged 
counterions that neutralize the charge of the macroions 
and 2ns salt ions, consisting half-and-half of counterions 
and oppositely charged colons. The inverse of the screen- 
ing parameter, the so-called Debye length i?D = 
"measures" the thickness of the neutralizing counterion 
layer around the macroions. Equation ^ shows that Rd 
can be varied by changing the properties of the solvent, 
in particular the salt concentration. 

Alexander et al. have demonstrated that many 
charged colloidal systems with highly charged macroions 
can be described to some extent by a Yukawa potential of 
the form of Eq. ^ , although the DH limit is restricted 
to particles with small charge. This is due to the fact 
that highly charged colloids have a strong tendency to 
form ordered structures at relatively low densities, i.e. at 
densities where the mean distance between neighboring 
macroions is much larger than their size. In such systems, 
each macroion has a very similar environment of mi- 
croions, and thus a reasonable approximation is to reduce 
the problem of computing the effective many-particle in- 
teractions between macroions to that of determining a 
mean-field potential of one particle in its Wigner-Seitz 
(WS) cell surrounded by a reservoir of salt ions Q. For 
spherical macroions, the WS cell is approximated by a 
sphere of radius R. Then, one considers the nonlinear PB 
equation for the single particle with appropriate bound- 
ary conditions , 



eps 



u[r) = 



Aire 



exp 



eu{r) 



exp 



eu{r) 



(3) 



n ■ Vw(r) 



Ze 



n ■ V-it(r) =0 {r = R) 



ksT 
{a/2 <r<R) 

(r = a/2) (4) 

(5) 



In this section, we consider charged spherical 
macroions of diameter a and positive charge Ze (here. 



with n the normal vector pointing outwards the sphere's 
surface and ps the salt concentration in the reservoir. 
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Equations and jSJ are solved numerically. Then, one 
assumes that at the cell boundary, i.e. far away from 
the surface of the particle, the solution u{r) of the PB 
equation can be approximated by an effective Yukawa 
potential. 



Zcff exp(Kcffg/2) 

1 + KcffO'/2 



exp(-Kcffr) 



(6) 



The parameters Zcs and KcS can be fixed by matching 
the effective potential at the cell boundary (r — R) 
with that of the solution of the nonlinear PB equa- 
tion. In the original paper by Alexander et al. @, 
this is achieved by the following recipe: The screen- 
ing parameter is determined by the microion den- 
sities n!l — po exp (±^) at the WS cell boundary via 

K^g — 47rAB(7T.:^ -f- n^). The effective charge Z^b is fixed 
as follows: First, Eqs. and ^ are linearized at r = i?. 
For the linearized equation a solution is determined such 
that the linear and the nonlinear solution match up to 
the second derivative at the cell boundary. Finally, Z^s 
is calculated from the integral over the charge density 
associated with the linear solution. 

A more elegant recipe to obtain Z^^ and HcS has re- 
cently been proposed by Trizac et al. Q . They show that 
the full numerical solution of the nonlinear PB equation 
is not needed to estimate the latter parameters. Instead, 
only the solution m/j at the cell boundary is required. 
Thus, only linearized equations have to be solved, and 
this can be done analytically. For Z^g, Trizac et al. 
find 

Z,g = ^^[{KlgaR/2 - 1) sinh(«:eff(i? - cj/2) (7) 
+ Koff(i?" cr/2)cosh(Koff(i?- cr/2)], 

where 70 = tanh(iii^). Equation implies Zcg < Z, 
where the effective charge (also called "renormalized 
charge" ) can be an order of magnitude smaller than the 
bare charge. 

The effective screening parameter is related to the 
effective salt concentration and Z^g via 



nf/V 



SttAb 



(l-7o)(l-'?)-^^cff"c(l-7o), (8) 



where represents the number of macroions per WS 
cell. The physical interpretation of the latter equation is 
related to the so-called Donnan effect. Since a macroion 
occupies a finite volume inside its WS cell, the microions 
are partially expelled. Thus, the salt concentration out- 
side the WS cell can be higher than inside the colloid 
compartement, or, in other words, there is an effective 
salt concentration which is smaller than the actual one, 
i.e. rig*^ < Ug. 

For dilute systems, where the Donnan effect should not 
be relevant, i.e. for 70,77 0, the effective salt concen- 
tration matches the actual one, and hence Eq. (jS)) can be 
rewritten as 



K-2 



47rAB(^cffnc + 2ns)/F. 



(9) 



Thus, setting = Ug, leads back to an one-parameter 
problem. It follows from Eq. that for monovalent 
microions, Z^g — Z also implies KcS = 

We emphasize that the systems considered in the fol- 
lowing do not match with the assumptions made in 
Alexander's concept of charge renormalization. In this 
work, "non-bulklike" systems are considered, for which 
the definition of a WS cell is not meaningful. Moreover, 
whereas charge renormalization is applied to distances far 
away from the surface of a macroion, we are interested 
in relatively small distances between macroions and thus 
also in the potential of mean force close to their surfaces. 



III. DETAILS OF THE SIMULATION 

Using classical MD simulations, we study charged col- 
loidal suspensions in the framework of the so-called prim- 
itive model. We consider systems of two or three posi- 
tively charged macroions of valency Z = Z^ = 255 and 
monovalent microions of charge Zcte — —1 (counterions) 
and of charge ZcoE = +1 (colons). The interaction po- 
tential between an ion of type a and an ion of type /3 
(a,/3 = m, ct, co), separated by a distance r from each 
other, is given by 

Z Z e^ 

Ual3 = " ^ h Aap exp{-Baf3{r ~ aap)l(Taf3} , (10) 

47rer 

where the dielectric constant is set to e = 79eo (with 
eo the vacuum dielectric constant), which corresponds 
to the value for water at room temperature. The pa- 
rameter (Tq^ is the distance between two ions at contact, 
— Ra + Rf3, where Ra is the radius of an ion of type 
a. In our simulations, we used i?ni = 10 nm and i?ct = 
Rco ~ 0.01i?in- The choice of the latter values guarantees 
that depletion effects are not relevant. The exponential 
in Eq. H1U|) is an approximation to a hard sphere inter- 
action for two ions at contact. For the parameters Aap 
we chose A^^ 1.84 eV, A^^t = ^mco = 0.0556544 eV, 
and Actct = ^ctco = ^coco = 0.0051 eV. The param- 
eters i?Q,/3 are all set to 3. The long-ranged Coulomb 
part of the potential and the forces were computed by 
Ewald sums in which we chose a — 0.05 for the constant 
and a cutoff wavenumber kc — 27r-\/66/ L in the Fourier 
part [l7j |. The linear dimension L of the simulation box 
is L = 159.026 nm, using periodic boundary conditions. 

Since the potential, Eq. ((117)) . is long-ranged, one has 
to consider the possible emergence of finite-size effects. 
However, in our simulations, the distance of a macroion 
to its next periodic image was always larger than 7 a 
(with a = tTmm), and, as discussed in the next section, 
at this distance the Coulomb interaction is sufficiently 
screened. Instead of using periodic boundary conditions, 
an alternative approach would be to confine the system 
by walls [HEl. However, due to the interaction of the 
ions with the walls, also in this case finite-size effects are 
relevant (indeed in Refs. [Tll ll^ . a correction term had 



FIG. 1: Electric field around macroion as a function of dis- 
tance for indicated salt concentrations. Solid lines are fits to 
Eq. @, where Zcft and Kcff are used as fit parameters. Data 
is shifted such that the dotted lines represent an electric field 
of -E = for each value of ng. Stars indicate critical macroion 
separations, as defined in text. Dashed lines are guides to 
the eye. Statistical errors are smaller than twice the size of 
symbols. 

to be introduced to estimate the "bulk" effective force 
between macroions). 

In order to determine the effective forces between 
macroions at a distance r, the macroions are fixed by 
decoupling macroionic and microionic time scales. This 
is achieved by assigning a mass to the macroion which is 
fO^mct (with TTict = 'Tico the mass of the microions). All 
the simulations were done at the temperature T = 298 K. 
Thus, the Bjerrum length for our system is Ab ~ 0.7Inm. 
The number of added colons was varied from rig = to 
rig = 1280. Being n'^ the number of macroions in the 
system, charge neutrality requires ritot = Zn'^ + 2ns for 
the total number of microions. 

The equations of motion were integrated using the ve- 
locity form of the Verlet algorithm. The simulations were 
done at constant temperature. In order to thermostat the 
system, it was coupled to a stochastic heat bath 17]. For 
a given set of parameters (rig, macroion separation r), we 
examined at least three independent start configurations. 
For equilibration, runs of 10^ — 10^ MD time steps were 
done, followed by a similar number of steps to calculate 
time averages. Depending on salt concentration, the time 
step varied from Ai = 1 • 10~''to to A< = 3 • 10~^to [with 

To = Rct^Jmct/{kBT)]. 



IV. RESULTS AND DISCUSSION 

A. Two macroions 

The effective interaction between two macroions can be 
quantified by the electric field E{r) around a macroion, 
depending on the separation r between the macroion's 
centers. The field E is given by E{r) = -^F{r), where 

F{r) = —^^7"^ is the total force on a macroion projected 
onto the line that connects the macroion's centers. Fig- 



FIG. 2: Logarithm of electric field times macroion separation 
versus r for various salt concentrations. Solid lines are fits to 
Eq. © using the same fitting parameters as in Fig. Q Data 
for ris = is shifted by 4, ria = 160 by 3, Ua = 320 by 2, and 
ris = 640 by 1, respectively. Stars and dashed lines have same 
meaning as in Fig. 



ureHshows E{r) for different salt concentrations, as indi- 
cated. A comparison to the prediction from DH theory, 
using Ucs{r) from Eq. with Z^ff and Kcff as fitting 
parameters, reveals a good agreement. However, as we 
will show in the following, the numerical values we ob- 
tain for Zcs and Kcff can neither be described by the DH 
predictions nor by the charge renormalization concept, 
as it results in Eqs. iQ and ||SJ). Plotting our data on a 
logarithmic scale exposes deviations from the DH form 
(Fig. |2Jl, which seem to increase with distance and salt 
concentration. We confirmed the absence of finite size 
effects by repeating some of our simulation runs in a sys- 
tem of double volume, finding our results to be consistent 
and not depending on the volume. 

To analyze our data further, it is useful to quantify the 
extent to what the Debye layers around the macroions 
overlap. To this end, we use the DH expression, Eq. jSJ, 
to estimate the inverse Debye layer thickness. According 
to this definition of k, Debye layers overlap, if -^(^ — 1) < 

1. For a given value of rig we thus define ri]) = a+2K^^ as 
an upper critical macroion separation. ForKcr(^ — 1) < 1, 
the macroions are (partially) located within each other's 
Debye layers. For a fixed salt concentration, we therefore 
define = a + k^^. Finally, we consider the limit, 
where a macroion's center of mass is located within the 
Debye layer of the other macroion. This is the case for 
Kcr(^ — ^) < 1, or. Tot = + The radu Tov^''^'' are 
indicated in Figs. ^ and 13 as stars (connected by dashed 
lines as a guide to the eye). 

As we see in Figs. ^ and |2 Debye layers overlap for 
almost all parameter combinations considered. Thus, 
we find a Yukawa-like potential also in the region of 
strongly overlapping Debye layers. The Yukawa-form 

(3) 

persists even for macroion-separations r < rov , where the 
macroions are relatively close to contact. These findings 
are in agreement with numerical solutions of the nonlin- 
ear PB equation for a system of two macroions . The 
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FIG. 3: a) Effective charge ZcS normalized by the bare charge 
Z and b) effective screening parameter KcB divided by k from 
DH theory [Eq. as a function of the number of salt ion 
pairs ris. The parameters ZcB and k^h result from fits shown 
in Fig.0and Fig.|21that include, as indicated, all data points 
with r > 1.3(T or only those with r > 1.82a. 
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FIG. 4: Test of Eq.@ in the limit 77,70 0. 



significant deviations from the DH fits seem to occur for 
non-overlapping Debye layers (see Fig.|21). This might be 
related to the low signal-to-noise ratio, which becomes 
worse for increasing values of Kr. 

Fig-Elshows the fitted values of Z^g and func- 
tion of salt ion pairs Ug- Here, we have normalized Z^g by 
the bare charge of the macroions and Keff by the screen- 
ing parameter k as predicted by Eq. The fit values 
in Fig. O are extracted from two different types of fits. In 
addition to fits that include all available data points, we 
performed also fits that are restricted to data points with 

(3) 

macroion separations of r > 1.82a > rov • Thus, in the 
latter fits, we exclude distances for which the center of 
a macroion penetrates into the Debye layer of the other 
one. 

As one can infer from Fig. |31 / k and Zcs /Z de- 
viate significantly from unity for all salt concentrations 
considered (except for Us ~ 250 where Kcs/k is close to 
one). Thus, for the systems under consideration, DH 
theory does not correctly describe the effective interac- 
tions. This indicates that a combination of nonlinear 
effects and, possibly, microionic correlations are relevant 
in the present case. 

Of particular interest is the behavior of ZcS- For all 
salt concentrations it is larger than the bare charge and 
it tends to increase with increasing salt content. This 
effect is even more pronounced if the fits are restricted 
to macroion separations of r > 1.82ct. This finding is in 
disagreement with the concept of charge renormalization, 
where one expects a decrease of the effective charge with 
increasing salt concentration. Indeed, cell model calcu- 
lations using the parameters of our MD simulations lead 
to Zeff « 0.8Z I13. 

A failure of the charge renormalization concept in the 
systems considered here is not surprising. For an isolated 
pair of macroions, there is no meaningful definition of a 
WS cell. Hence, it is difficult to define the volume frac- 
tion 7] reasonably. Provided, that the WS cell can be 
considered as a sphere around a macroion, we can come 
up with two boundaries: The WS cell should not pene- 



trate the Debye layer, thus, R > ri]) /2 = + a/2. In 
addition, the two WS cells should not overlap, hence, 
R < r/2. Note, that the upper boundary depends 
on the macroion separation, whereas, according to the- 
ory, K is not a function of r. However, our simula- 
tion data indicates, that -^KcS is slightly different from 
zero. Taking into account both limits of i?, the vol- 
ume fraction rj = "'c(^)'^ should fullfill the inequality 
(1 -f < < f. With K from Eq. Q, the lower 

boundary takes values from 1.6-10"^ (rig = 0) to 9.3-10"^ 
{ug — 1280). The overall volume fraction of macroions 
in our simulation box is given by rj' — ~ 2.I • 10~^. 
Hence, in dilute systems, rj' cannot be regarded as the 
relevant volume fraction for testing Alexander's charge 
renormalization concept. Moreover, 77' ^ 1 indicates, 
that Eq. Q should hold at least for small salt concentra- 
tions, provided that 70 can be neglected, i.e., Kcff ~ 
Although both values, Zcs and KcS, are not found to 
be in agreement with DH theory, their relation seems to 
be compatible with Eq. @, as can be seen from Fig. ^ 
Thus, to a good approximation, the effective salt concen- 
tration matches the actual one, i.e., the Donnan effect is 
indeed negligible. 

So far, we have addressed only the behavior of effec- 
tive pair forces. In order to analyze the microionic de- 
grees of freedom, we calculate the angular resolved nega- 
tive charge density distribution, p_(a), which we define 
as follows: For a given macroion separation, we draw a 
sphere of radius r/2 around each macroion and project 
all counterions, which are located within this fictitious 
"WS cell", onto a plane, which contains the macroions' 
centers. Being nct{a) the number of counterions at angle 
a, which is taken relative to the connecting line to the 
other macroion, p-{a) is given by 

p_(a) = 10-Vt(a)(^)'. (11) 

Note, that P-{a) is normalized via the volume of the WS 
sphere, nr^/G. 
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FIG. 5: Angular resolved negative charge density distribution 
around macroion, as defined by Eq. (IIH . for the indicated 
salt concentrations and macroion separations. Data sets are 
shifted by 0, 0.05, 0.1, and 0.15 (from below). The solid lines 
are calculated from the superposition of one-particle charge 
density distributions from DH theory (for details see text). 




FIG. 6: Deviation from charge neutrality within "WS cell", 
measured via Q{R), as defined in Eq. 1131 . Data is plotted for 
Tia — 1280. Results are compared to Eq. 1)16^ . using k from 
Eq. J^J (solid line) and Atcff (dashed line). 



Fig. [SI shows P-{a) for various combinations of nr. 
Four different cases are considered: For rig — 0, we choose 

For Us — 320, we consider the 
For ris = 1280, we take 



(3) 



r such that Hcsr <SC nr, 

(3) (2) 

case HiToY ^ ^cfi^ ^ /^^ov • 



Kr^ < Kcsr < /trov"* and Hesr ^ Kriv\ respectively. 

For an isolated pair of macroions, the electric field 
around a macroion only exhibits a spherical symmetry 
in the limit Kr — > oo. Therefore, one might expect that 
P-{a) is not independent of a, and thus, it should reveal 
anisotropics. However, within the statistics of our data, 
we find flat distributions within the "WS cell" . This 
holds even for the smallest value of nr, where the Debye 
layer around a given macroion is strongly perturbed by 
the other macroion. The occurrence of isotropic distri- 
butions P-(ck) might be due to nonlinearities, which are 
of course not accounted for in the DH limit. In order 
to rationalize this hypothesis, we checked whether p-{a) 
can be "reconstructed" by a naive superposition of coun- 
terion charge distributions around a single macroion. To 
this end, we considered first such single-particle distri- 
butions as obtained from DH theory using the screening 
parameter k as given by Eq. Q and the bare charge Z 
for the charge of the macroion [note that the charge dis- 
tribution is just proportional to the potential given by 
Eq. © in the linearized DH limit]. Then, we projected 
the superposition of the latter distributions onto a cubic 
lattice with 10* grid points. From this, we finally calcu- 
lated p-{a). The results are included in Fig. [S] as solid 
lines. We clearly see that the so calculated P-{a) are 
anisotropic, and this indicates that the flat distributions 
obtained from the MD simulations might be due to the 
occurrence of nonlinearities. 

The behavior of P-{a) might also explain why the ef- 
fective charge Z^s is higher than the bare charge Z, in 
contrast to the prediction from the charge renormaliza- 
tion concept. We can infer from the angular distribu- 



tions p-{a) that there are less counterions between the 
macroions than expected from a naive superposition prin- 
ciple. This effect might be of entropic origin indicating 
that the entropy gain related to isotropic distributions 
dominates over energetic contributions. However, ener- 
getically unfavoured microion distributions might yield 
an additional repulsion between the macroions, and this 
might explain the flnding that Z^g is larger than the bare 
charge. 

We have already mentioned that the introduction of 
a "WS cell" is not appropriate for a system of two iso- 
lated macroions and thus cell models that lead to charge 
renormalization cannot be applied. There is also another 
reason why the concept of charge renormalization is not 
appropriate in the present case. If we define the bound- 
ary of the (spherical) "WS cell" by the sphere of radius 
R — r/2 around a macroion, this cell is not a neutral ob- 
ject, i.e. the total charge inside the cell is nonzero. This is 
in contrast to the assumptions of Alexander's cell model 
which is not applicable for small macroion separations. 

However, it is instructive to study the total charge of 
the "WS cell" for our system. Charge neutrality requires 



\2R) 



dap^{a) = / dap_(a). 

"'0 



(12) 



where p+{a) has an analogous definition as p_(a), but 
now the number of counterions at angle a is replaced 
by the corresponding number of colons. It follows from 
Eq. lO that 



Q{R) = 1 



1 f2R 



Z \ G 



da[p+(a)-p_(a)] (13) 



should vanish, if the overall charge within the "WS cell" 
is zero. In Fig. we show (5(i?) for a fixed salt con- 
centration of rig = 1280. It is compared to an estimate, 
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which follows from DH theory: Suppose, the counterion 
density around a macroion is given by an expression of 
the DH form, 



p(r) 



7exp[~K(r - f ) 
(1 + if )rA| 



(14) 



where 7 is a dimensionless normalization constant. Since 
Q{R) = ^ ~ Iv{r) d'^^p(r), the total charge reads 



QiR) = l-^ r drr^p{r) 



= 1 - 
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(15) 



The normalization constant introduced in Eq. (|14(l is de- 
termined by the boundary condition Q{R 00) — > 0, 

thus, 7 = ^'^^/2k+i/k'^) ■ Hence, the charge inside the 
fictitious WS cell of radius R is given by the expression 



QiR) = 



1 + kR 
^ ^ 2 



exp 



r fR 




— Ka[ 




V (7 





(16) 



Note, that the second boundary condition, Q{R= ^) = 
1, is intrinsically fulfilled. If we identify the inverse 
screening length as k from Eq. Eq. (|16|) slightly un- 
derestimates Q{R). Replacing k, by the effective inverse 
screening length leads to an almost perfect agreement 
with the results of our simulations. 
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FIG. 7: Relative deviation A of electric field around outer- 
most macroion in the coaxial geometry as compared to expec- 
tation for pairwise additivity for different distances ri2 and 
r23 between the particles. 



Here, i?*^^) is the superposition of the two-body inter- 
actions calculated in the previous subsection, whereas 
E^^'^ = —-^■§pV{r) follows from the force acting on the 
outermost particle in the three-macroion configuration. 

For the coaxial geometry, the pair contribution is given 
by E^"^^ — —^-^ J2i>i ^(''li)- the configuration of an 
equilateral triangle with side length R, one has to take 
into account, that the forces do not act along the same 
direction. If we denote the positions of the macroions by 
Ri [i — 1,2,3) the effective force F{R) is given by the 
total force Fi on particle 1 projected onto the difference 
vector d= Ri -\{Ri+R2+ R3), 



B. Three macroions 

In this section, we consider systems with three 
macroions in two different geometries by placing them 
along a line or at the corners of an equilateral triangle. 

In general, the interaction energy for three particles 
can be written as 

V{r) = T/i2(ri2)-l-Vi3(ri3)-f V'23(r23) + Vi23(ri23), (17) 

where Vij (r.y ) is the pair potential between particle i and 
j. The last term on the right hand side represents the 
three-body interactions. 

We measure the force on the outermost particle (par- 
ticle "1") and define the relative deviation of the electric 
field with respect to the expectation for pairwise additiv- 
ity by 

^ i:[V{r)-V{n2)~V{n3)] 



A = - 



FiR) = Fi • ^ 
\d\ 



(19) 



(18) 



Thus, the two-body contribution in the equilateral tri- 
angle is E(^^ = -||: E,>i Viru) cos(^/6). 

Results for the coaxial geometry are displayed in Fig.[7| 
where the deviation A from pairwise additivity is plotted 
as a function of the parameter /^^ = "lik [ ilialiia) _ 
1]. The quantity /^^ describes the overlap between the 
Debye layers around the three macroions. For /^^ < 1, 
the three Debye layers exhibit an overlap. 

One can infer from Fig. [7|that the three-body interac- 
tion between the macroions yields attractive corrections 
to pairwise additivity. At small salt concentration, i.e. if 
/o^ is significantly smaller than one for a given distance 
between the macroions, three-body corrections are most 
pronounced and they are weakly dependent on /^^ . But 
if reaches values that are of the order of one, the pa- 
rameter A increases rapidly and seems to vanish at high 
values of f^^. Thus, three-body contributions are of im- 
portance if there is an overlap between the three Debye 
layers. This shows that the range of three-body contri- 
butions is of the order of the Debye length and thus the 
concept of screening is also very useful for the discussion 
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FIG. 8: Relative deviation of electric field around macroion as 
compared to expectation for pairwise additivity versus overlap 
factor /ov (see text). Triangular symbols represent triangular 
setup, circles represent coaxial geometry. 
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FIG. 9: Negative charge density distribution for three- 
macroion case, comparing the coaxial geometry to the tri- 
angular configuration. The latter is shifted by —0.05. The 
macroion separation is fixed to r = 1.3a. Salt concentration 
is Ws = (left) and ris — 1280 (right), respectively. Solid lines 
have same meaning as in Fig. |3 



of many-body effects. 

It is interesting that the three-body terms are much 
smaller in the coaxial geometry if the distance between 
neighboring macroions is close to contact. This corre- 
sponds to the data for ri2 = r23 = 1.3ct in Fig. 13 In this 
case, the pair interaction probably is the most important 
contribution because the three-body force on the outer- 
most particle is effectively screened out by the particle in 
the middle. 

In the case of an equilateral triangle, the condition for 
three overlapping Debye layers reads /^^ = k(^V(-^^ — 

i) < 1. Using the same values of k and r for both ge- 
ometries leads to < /q^, thus, in the triangular con- 
figuration, the Debye layers exhibit a stronger overlap. 

In Fig. |S1 three different geometries for the macroion 
triple are compared: 

1. an equilateral triangle with side length r — 1.3cr 



2. an equilateral triangle with side length r = 1.82a 

3. the coaxial geometry with r 



ri2 = r23 



1.3a. 



The strongest triplet interactions are revealed for case 
1. Different from the coaxial geometry, the magnitude 
of the parameter A increases with decreasing distance r 
between the particles in the triangle. This can be eas- 
ily understood since the interaction between the three 
macroions in the triangular geometry is not effectively 
screened by one of them but only by the microions in the 
Debye layers. For case 1 the magnitude of the parameter 
A also seems to increase with increasing /*^. Up to now, 
we do not have an explanation for this behavior. 

Comparing case 2 and 3, we see that deviations from 
pairwise additivity are similar in both cases, and the over- 
lap parameters and /^^ have comparable values. A 
similar feature has been found in a numerical solution 
of the nonlinear PB equation by Russ et al. 0. These 
authors report that the three-body potential is indepen- 
dent of geometry if the sum over the distances between 
neighboring particles is constant. 

We would like to point out that there is always a triv- 
ial contribution to the many-body potential, which stems 
from an increased microion concentration associated with 
the addition of macroions. Thus, for a fixed volume, 
the effective screening length of the system is decreased. 
A comparison between the measured three-body term, 
E'^'^\ and the two-body contribution, E'^^\ taken from 
the pure two-macroion case should therefore in general 
yield a non-pairwise additivity. From that point of view, 
DH theory already predicts a correction to pairwise ad- 
ditivity. 

Similar to the previous subsection, we calculate the 
angular resolved charge density distribution around the 
outermost macroion [see Eq. In order to account 

for the differences between coaxial and triangular geom- 
etry, we introduce an angle ao, such that the system is 
symmetric around a = ao. Thus, we have to choose 
ao — for the coaxial geometry and ao = 7r/6 for the 
equilateral triangle. 

As Fig. |51 shows, the charge distribution for the three 
macroion case is not isotropic. Compared to the case 
of two macroions, the internal energy of the system is 
now larger, and, therefore, entropic contributions might 
be less important. In accordance with our previous con- 
sideration, namely that an energetically unfavored mi- 
croion distribution might lead to an additional repulsion 
between the macroions, we might conclude here that the 
onset of anisotropy is correlated with non-pairwise ad- 
ditivity. The choice of the same value of r leads to a 
stronger anisotropy for the triangle as compared to the 
coaxial geometry, which is consistent with the behavior 
of the parameter A (see Fig. |Sl . As done in the previ- 
ous section, we compare our charge distribution to the 
one which follows from naive superposition of the DH 
distributions. Surprisingly, for the three macroion con- 
figurations, this superposition seems to work very well. 
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V. CONCLUDING REMARKS 

We performed classical MD simulations in order to 
investigate effective interactions between isolated pairs 
and triples of charged macroions in the framework of the 
primitive model. 

On the pair level, these interactions are surprisingly 
well described by the DH limit of the PB equation. In 
particular there is no evidence for charge renormaliza- 
tion as predicted by cell models. These models would 
predict an effective charge which is considerably smaller 
than the bare charge of the macroion 0, 0, E3- This 
finding is not due to finite size effects in the simulation 
which might emerge if the Debye length exceeds the size 
of the simulation box. It rather follows from the fact, 
that the cell model must not be applied to systems of iso- 
lated macroions. This means that the concept of charge 
renormalization might be relevant for bulk systems, but, 
in the case of systems of isolated macroions, simulations 
should be compared to direct solutions of the nonlinear 
PB equation. 

In this work, we have studied systems with small salt 
concentrations of the order of a few /xMol, and we have 
considered configurations for which the Debye layers of 
the different macroions exhibit a strong overlap. An in- 
teresting result of our simulations is the occurrence of 
repulsive corrections to DH theory: For the macroion 
pair, we find effective macroion charges that are slightly 
higher than their bare charge. Similar results have been 
reported in previous work, e.g ., in an ab initio density 
functional theory approach , where the ratio between 
effective and bare charge was found to be between 1.06 
and 1.38, depending on salt concentration and the value 
of the bare charge. An "repulsive correction to DH the- 
ory" is also indicated by the isotropic charge distribution 
around the macroion in the case of the macroion pair. 
Such an isotropic distribution is not expected from a 
naive superposition of one-particle density distributions 
as obtained from DH theory. Hence, there seem to be less 
counterions between the macroions than expected from 



DH theory which can be related to an increase of the 
effective charge. The microscopic origin of this effect is 
not clear, but it might be of entropic origin. 

In agreement with previous analytical [loj . numerical 
^Ji, JJs J^J and experimental studies of systems with 
three isolated macroions, we find that corrections to non- 
pairwise additivity (and thus the three-body terms in the 
effective potential) are attractive. The strength of these 
attractive contributions is strongly correlated with the 
overlap of all three Debye layers. This shows that the 
concept of a screening length is also very useful to quan- 
tify the effect of three-body interactions. Different from 
the case of two macroions, the charge distribution in the 
three-macroion case is anisotropic. In this case, the sim- 
ple superposition of three one-particle density distribu- 
tions from DH theory yields a rather good description of 
the charge distribution in the three-macroion case. This 
finding seems to agree with a recent numerical solution 
of the PB equation for three isolated macroions [l6l |. 

In further simulation studies, we will investigate inter- 
actions between more than three particles to understand 
the crossover to bulk effective interactions. In the lat- 
ter case, the concept of charge renormalization seems to 
work very well. Our present simulations suggest that 
many-body interactions in bulk systems yield renormal- 
ized charges that can be much smaller than the bare 
charges of charged colloidal particles. A profound under- 
standing of these issues might also pro vide new insight 
into electrophoresis experiments jl9ll20j |. 
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